

oo::class create Traj {
    constructor { Operator {filelist {} }} {
	my variable operator
	set operator $Operator
	my read $filelist
    }
    method eachframe { block } {}
    method every { timeperioid block} {}
    method save fname {}
    method read filelist {}
}

source trajAna.tcl
#core utilities
oo::define Traj method eachframe { block } {
    my variable operator
    my variable dcdlist
    foreach dcdfile $dcdlist {
	set path [$operator search $dcdfile]
	if { [string length $path] <= 0 } { error "$dcdfile not exists!" }
	$operator sendCmd "open read unit 31 file name $path"
	set result [$operator sendCmd "traj query unit 31"]
	set nframe 0
	puts "opened $path"
	regexp -all {NUMBER\s+OF\s+COORDINATE\s+SETS\s+IN\s+FILE:\s+(\d+)\n} $result  result nframe
	$operator sendCmd "trajectory firstu 31 nunit 1"
	for {set i 0} { $i < $nframe } { incr i}  {
	    $operator sendCmd "traj read"
	    uplevel 1 "$operator do {$block}"
	} 
	$operator sendCmd "close unit 31" 
    }
}

oo::define Traj method every { timeperoid block} {
    my variable operator
    my variable dcdlist
    #query the dcd file for 
    set nskip 1
    foreach dcdfile $dcdlist {
	set path [$operator search $dcdfile]
	if { [string length $path] <= 0 } { error "$dcdfile not exists!" }
	$operator sendCmd "open read unit 31 file name $path"
	set result [$operator sendCmd "traj query unit 31"]
	set nframe 0
	puts "opened $path"

	regexp -all {NUMBER\s+OF\s+COORDINATE\s+SETS\s+IN\s+FILE:\s+(\d+)\n} $result  matched nframe
	regexp -all {FREQUENCY\s+FOR\s+SAVING\s+COORDINATES:\s+(\d+)\n} $result matched savefrq
	regexp -all {THE\s+INTEGRATION\s+TIME\s+STEP\s+\(PS\):\s+(\d+\.\d+)\n} $result matched timesteps
        set interval [expr $timesteps *$savefrq]
        #puts "interval $interval"
        switch -glob $timeperoid {
          *ps { if {[scan $timeperoid "%dps" t ] >= 1} {set nskip [expr int($t/$interval)]}; }
          *ns {if {[scan $timeperoid "%dns" t ] >= 1} {set nskip [expr int($t/$interval)]};}
          *step {if {[scan $timeperoid "%dstep" t ] >= 1} {set nskip [expr int($t)]};}
         default { set nskip [expr int($timeperoid)]}
         }
	$operator sendCmd "trajectory firstu 31 nunit 1"
	for {set i 0} { $i < $nframe } { incr i }  {
	    $operator sendCmd "traj read"
	    if {[expr $i % $nskip] == 0 } {
	    	uplevel 1 "$operator do {$block}"
	    }
	} 
	$operator sendCmd "close unit 31" 
    }
}

oo::define Traj method save fname {
     my variable operator
     my variable dcdlist
     set nfiles [llength $dcdlist]
     set stepsize 50
     for {set i 0 } { $i < $nfiles } {incr i $stepsize} {
       set flist [lrange $dcdlist $i [expr $i +$stepsize]]
       if {$i == 0 } {
         $operator sendCmd "coor copy comp"
         my merge2 $flist tmp.dcd
       } else {
         my merge2 [concat tmp.dcd $flist] tmp2.dcd
         exec mv tmp2.dcd tmp.dcd 
       }
     }
     exec mv tmp.dcd $fname
}
#setup 
oo::define Traj method read filelist {
    my variable dcdlist
    set dcdlist $filelist
}
    

#Helper functions
oo::define Traj method query dcdfile {
    my variable operator
    $operator sendCmd 	""
}
# maximun merge 50 files from discontinous simulation
oo::define Traj method merge2 { flist fname } {
    my variable operator
    set nfile [llength $flist]
    $operator sendCmd "bomblev -2"
    for {set i 0} {$i <$nfile} {incr i} {
       	set path [$operator search [lindex $flist $i]]
	if { [string length $path] <= 0 } { error "[lindex $flist $i] not exists!" }
	$operator sendCmd "open read unit [expr 30+$i] file name $path"
    }
    $operator sendCmd "open write unit 29 file name $fname"
    $operator sendCmd "trajectory firstu 30 nunit $nfile iwrite 29 "
    for {set i 0} {$i <$nfile} {incr i} {
	set result [$operator sendCmd "traj query unit [expr 30+$i] "]
	set nframe 0
	regexp -all {NUMBER\s+OF\s+COORDINATE\s+SETS\s+IN\s+FILE:\s+(\d+)\n} $result  result nframe
	for {set j 0} { $j < $nframe } { incr j}  {
	    $operator sendCmd "traj read"
	    # centering and orienting the peptide
	    $operator sendCmd "coor orient rms" 
            $operator sendCmd "traj write"
	} 
    }
    
    for {set i 0} {$i <$nfile} {incr i} {
       $operator sendCmd "close unit [expr 30+$i]"
    }
    $operator sendCmd "close unit 29"
    $operator sendCmd "bomblev -2"
}

# maximun merge 9 files from continous simulation
oo::define Traj method merge { flist fname {selection "select all end"}  } {
    my variable operator
    set nfile [llength $flist]
    for {set i 0} {$i <$nfile} {incr i} {
       	set path [$operator search [lindex $flist $i]]
	if { [string length $path] <= 0 } { error "$dcdfile not exists!" }
	$operator sendCmd "open read unit [expr $i+30] file name $path"
    }
    $operator sendCmd "open write unit 29 file name $fname"
    $operator sendCmd "merge firstu 30 nunit $nfile outp 29 skip 1 -\n $selection "
    for {set i 0} {$i <$nfile} {incr i} {
       $operator sendCmd "close unit [expr 30+$i]"
    }
    $operator sendCmd "close unit 29"
}

# utilities
proc histogram { l start end nbins } {
    foreach v $l {
	set i [ expr floor(($v-$start)/$nbins)]
	if [info exists hist($i)] {
	    set hist($i) [expr $hist($i) +1]
	} else {
	    set hist($i) 1
	}
    }
    return [array get hist]
}
proc groupby { ic filters } {
    upvar $ic IC
    for {set i 0 } { $i < [llength $IC(phi)]} { incr i} {
	set group ""
	foreach f $filters {
	    switch $f {
		residue {
		    set atom [lindex $IC(resid) $i]
		    append group [belongtoResidue $atom]
		    append group ,
		}
		backbone {
		    set atom [lindex $IC(atom) $i ]
		    append group [belongtoBB $atom]
		    append group ,
		}
	    }
	}
	set group [string range $group 0 end-1]
	if [info exists Group($group)] {
	    lappend Group($group) [lindex $IC(num) $i]
	} else {
	    set Group($group) [list [lindex $IC(num) $i]]
	}
    }
    return [array get Group]
}
proc belongtoBB { atoms } {
    set BackboneAtom {CA N C O *C *CA *N *O }
    if [contain_only $atoms $BackboneAtom ] { return BB 
    } else {
	return SC
    }
}
proc belongtoResidue { atoms } {
    foreach atom $atoms {
	if [info exists num($atom)] {
	    set num($atom) [expr $num($atom) +1]
	} else {
	    set num($atom) 1
	}
    }
    foreach atom [array names num] {
	if [info exists maxatom] {
	    if { $num($atom) >= $num($maxatom) } {
		set maxatom $atom
	    }
	} else {
	    set maxatom $atom
	}
    }
    return $maxatom
}
proc contain_only { table Set } {
    set only 1
    foreach item $table {
	if { [lsearch $Set $item] <0} {
	    set only 0
	}
    }
    return $only
}
proc count { atomlist } {
    foreach atom $atomlist {
	if [info exists num($atom)] {
	    set num($atom) [expr $num($atom) +1]
	} else {
	    set num($atom) 1
	}
    }
}
